# Read in ACCESS full dataset --------------------------------------------------
# Note: This treats the two waves as repeated cross-sections

pool <- read.dta13(here("Data", "poolweights_dec31.dta"),
                   convert.factors = T)

# Set up data dataset for analysis ---------------------------------------------

pool <- pool  %>% 
  mutate(wave = as.numeric(year == 2018),
         caste = m1_q25_caste,
         scst = as.numeric(caste == "Scheduled Caste" | 
                             caste == "Scheduled Tribe"),
         hhsize = m1_q27_no_adults + m1_q29_no_children,
         lnexp = log(m1_q32_month_expenditure + 1),
         bplaay = as.numeric(m1_q26_ration == "BPL" | m1_q26_ration == "Antyodaya"),
         useelec = as.numeric(m2_q68_elec == "Yes"),
         usegrid = as.numeric(m2_q55_grid == "Yes"),
         uselpg = as.numeric(m4_q103_lpg == "Yes"),
         duration_day = m2_q69_elec_hrs,
         reliability = m2_q71_elec_out_days,
         elec_duration_tier = elec_duration_tier,
         elec_reliability_tier = elec_reliability_tier,
         satisfied = as.numeric(m2_q77_electrified_satisfaction == "Satisfied")
         ) %>% 
  mutate_at(vars(matches("m2_q75")),
            .funs = ~(as.numeric(. > 0 & !is.na(.)))
            ) %>% 
  mutate(ownincbulb = as.numeric(m2_q75_1_no_Incandescent_bulb == 1),
         owncflbulb = as.numeric(m2_q75_2_no_cfl_bulb == 1),
         ownledlight = as.numeric(m2_q75_3_no_led_light == 1),
         owntublight = as.numeric(m2_q75_4_no_tube_light == 1),
         ownmobile = as.numeric(m2_q67_1_phone_no > 0 & !is.na(m2_q67_1_phone_no)),
         ownradio = as.numeric(m2_q75_9_no_radios == 1),
         ownfan = as.numeric(m2_q75_5_no_fans == 1),
         owntv = as.numeric(m2_q75_8_no_tvs == 1),
         owniron = as.numeric(m2_q75_6_no_irons == 1),
         ownfridge = as.numeric(m2_q75_7_no_fridges == 1),
         owncooler = as.numeric(m2_q75_10_no_cooler == 1),
         ownwshmchn = as.numeric(m2_q75_11_no_wash_machine == 1),
         ownelstove = as.numeric(m2_q75_12_no_elec_stove == 1),
         owninvert = as.numeric(m2_q75_13_no_inverter == 1),
         ownelpump = as.numeric(m2_q75_14_no_elec_pump == 1)
  ) %>% 
  mutate(
         srv_light = as.numeric(
           rowSums(select(., ownincbulb, owncflbulb, ownledlight, owntublight), na.rm = T) >= 1),
         srv_ict = as.numeric(
           rowSums(select(., ownmobile, ownradio), na.rm = T) >= 1),
         srv_ent = as.numeric(
           rowSums(select(., owntv), na.rm = T) >= 1),
         srv_spch = as.numeric(
           rowSums(select(.,ownfan, owncooler), na.rm = T) >= 1),
         srv_frdg = as.numeric(
           rowSums(select(.,ownfridge), na.rm = T) >= 1),
         srv_mech = as.numeric(
           rowSums(select(.,ownwshmchn, ownelpump), na.rm = T) >= 1),
         srv_thrm = as.numeric(
           rowSums(select(.,owniron), na.rm = T) >= 1),
         srv_cook = as.numeric(
           rowSums(select(.,ownelstove), na.rm = T) >= 1),
         srv_light_num = rowSums(select(., m2_q75_1_no_Incandescent_bulb, 
                                        m2_q75_2_no_cfl_bulb, m2_q75_3_no_led_light, 
                                        m2_q75_4_no_tube_light), na.rm = T),
         srv_ict_num = rowSums(select(., m2_q67_1_phone_no, m2_q75_9_no_radios), na.rm = T),
         srv_ent_num = rowSums(select(., m2_q75_8_no_tvs), na.rm = T),
         srv_spch_num = rowSums(select(.,m2_q75_5_no_fans, m2_q75_10_no_cooler), na.rm = T),
         srv_frdg_num = rowSums(select(.,m2_q75_7_no_fridges), na.rm = T),
         srv_mech_num = rowSums(select(.,m2_q75_11_no_wash_machine, m2_q75_14_no_elec_pump), na.rm = T),
         srv_thrm_num = rowSums(select(.,m2_q75_6_no_irons), na.rm = T),
         srv_cook_num = rowSums(select(.,m2_q75_12_no_elec_stove), na.rm = T)
         ) %>%
  mutate(srv_tot = rowSums(select(.,srv_light:srv_cook)),
         srv_tot_num = rowSums(select(.,srv_light_num:srv_cook_num))
         )

pool <- pool %>% 
  select(state = m1_q8_state, district = m1_q9_district, village = m1_q11_village_code,
         hh = finalhhid, wave, caste, scst, hhsize, lnexp,
         exp = m1_q32_month_expenditure, bplaay, education = m1_q23_edu,
         religion = m1_q24_religion, housetype = m1_q39_house_type,
         useelec, usegrid, uselpg, duration_day, reliability,
         ownincbulb:ownelpump, srv_light:srv_tot,
         elec_duration_tier, elec_reliability_tier,
         ceewtier = elec_overall_tier, satisfied, weights) %>%
  mutate(educ10th = as.numeric(education == "Up to 10th Standard" |
                                 education == "12th Standard or Diploma" |
                                 education == "Graduate and Above"),
         state = as.character(factor(state)),
         puccahouse = as.numeric(housetype == "Pucca"),
         unelectrified = ifelse(useelec == 1, 0, 1)) %>% 
  droplevels() %>% 
  arrange(state, village, hh, wave)

# Merge in village dataset key variables ---------------------------------------

villimport2014 <- read.dta13(here("Data", "vill2014.dta"),
                             nonint.factors = T, generate.factors = T)

villimport2018 <- read.dta13(here("Data", "vill2018.dta"),
                             nonint.factors = T, generate.factors = T)

vill2014 <- villimport2014 %>% 
  mutate(wave = 0,
         leaderscst = as.numeric(v_q13_caste == "Scheduled Caste" |
                                   v_q13_caste == "Scheduled Tribe"),
         villelec = as.numeric(v_q21_elec_grid == "Yes"),
         villmg = as.numeric(v_q22_microgrid == "Yes"),
         villlpg = as.numeric(v_q25_lpg == "Yes")) %>% 
  rowwise() %>% 
  mutate(villprisch = sum(v_q34_1_primary_elec,v_q34_1_primary_unelec),
         villprischelec = v_q34_1_primary_elec / villprisch,
         villmidsch = sum(v_q34_2_middle_elec,v_q34_2_middle_unelec),
         villmidschelec = v_q34_2_middle_elec / villmidsch,
         villangwdi = sum(v_q34_3_aanganwaadi_elec, v_q34_3_aanganwaadi_unelec),
         villangwdielec = v_q34_3_aanganwaadi_elec / villangwdi,
         villdispen = sum(v_q34_4_dispensary_elec, v_q34_4_dispensary_unelec),
         villdispenelec = v_q34_4_dispensary_elec / villdispen,
         villrelplc = sum(v_q34_5_religion_elec, v_q34_5_religion_unelec),
         villrelplcelec = v_q34_5_religion_elec / villrelplc,
         villcomspc = sum(villprisch, villmidsch, villangwdi, villdispen,
                          villrelplc, na.rm = T),
         villcomspcelec = sum(villprischelec, villmidschelec, villangwdielec,
                              villdispenelec, villrelplcelec, na.rm = T) / 
           villcomspc) %>% 
  ungroup() %>% 
  select(state = v_q1_state_name, division = v_q2_division_name,
         district = v_q3_district_name, village = v_q6_village_code,
         leadercaste = v_q13_caste, leaderscst, wave, leaderreligion = v_q12_religion,
         villelec, villmg, villlpg,villprisch,villprischelec,villmidsch,
         villmidschelec, villangwdi, villangwdielec, villdispen, villdispenelec,
         villrelplc, villrelplcelec, villcomspc, villcomspcelec)

vill2018 <- villimport2018 %>% 
  mutate(wave = 1,
         leaderscst = as.numeric(v_q13_caste == "Scheduled Caste" |
                                   v_q13_caste == "Scheduled Tribe"),
         villelec = as.numeric(v_q21_elec_grid == "Yes"),
         villmg = as.numeric(v_q22_microgrid == "Yes"),
         villlpg = as.numeric(v_q25_lpg == "Yes")) %>% 
  rowwise() %>% 
  mutate(villprisch = sum(v_q34_1_primary_elec,v_q34_1_primary_unelec),
         villprischelec = v_q34_1_primary_elec / villprisch,
         villmidsch = sum(v_q34_2_middle_elec,v_q34_2_middle_unelec),
         villmidschelec = v_q34_2_middle_elec / villmidsch,
         villangwdi = sum(v_q34_3_aanganwaadi_elec, v_q34_3_aanganwaadi_unelec),
         villangwdielec = v_q34_3_aanganwaadi_elec / villangwdi,
         villdispen = sum(v_q34_4_dispensary_elec, v_q34_4_dispensary_unelec),
         villdispenelec = v_q34_4_dispensary_elec / villdispen,
         villrelplc = sum(v_q34_5_religion_elec, v_q34_5_religion_unelec),
         villrelplcelec = v_q34_5_religion_elec / villrelplc,
         villcomspc = sum(villprisch, villmidsch, villangwdi, villdispen,
                          villrelplc, na.rm = T),
         villcomspcelec = sum(villprischelec, villmidschelec, villangwdielec,
                              villdispenelec, villrelplcelec, na.rm = T) / 
           villcomspc) %>% 
  ungroup() %>% 
  select(state = v_q1_state_name, division = v_q2_division_name,
         district = v_q3_district_name, village = v_q6_village_code,
         leadercaste = v_q13_caste, leaderscst, wave, leaderreligion = v_q12_religion,
         villelec, villmg, villlpg,villprisch,villprischelec,villmidsch,
         villmidschelec, villangwdi, villangwdielec, villdispen, villdispenelec,
         villrelplc, villrelplcelec, villcomspc, villcomspcelec)

vill <- rbind(vill2014,vill2018) %>% arrange(state, division, district, village, wave)

vill <- vill %>% 
  mutate_cond(district == "JAJPUR", district = "JAJAPUR") %>% 
  mutate_cond(district == "KANNUAJ", district = "KANNAUJ") %>% 
  mutate_cond(district == "SIDDHARTH NAGAR", district = "SIDDHARTHNAGAR") %>%
  mutate_cond(district == "BULANDSHAHR", district = "BULANDSHAHAR") %>% 
  mutate_cond(district == "BIJNOR", district = "BIJNOUR")

pool <- left_join(pool, vill, by = c("state","district","village","wave"))

rm(vill,vill2014,vill2018,villimport2014,villimport2018)

# Add in village coordinates ---------------------------------------------------

# p_load(raster, sf)
# 
# access <- raster(here("Data", "2015_accessibility_to_cities_v1.0.tif"))
# 
# villcoord <- read_xlsx(here("Data", "villcoord_dec27.xlsx")) %>% 
#   dplyr::select(village = Village_code, Longitude, Latitude) %>% 
#   st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)
# 
# villcoord <- villcoord %>%
#   mutate(timetocity = raster::extract(access, ., na.rm = T),
#          timetocity = as.numeric(timetocity))
# 
# saveRDS(villcoord, here("Data", "timetocity.rds"))

villcoord <- readRDS(here("Data", "timetocity.rds"))
 
pool <- left_join(pool, villcoord, by = c("village"))

rm(villcoord)
